### RELICATION SCRIPT for:                                            ###
### Thinking Ideologically: The Limited Role of Left and Right Labels ###
### as Policy Shortcuts                                               ###
### Authors: Lachance, Sarah & Clareta, Treger                        ###
### Journal: Public Opinion Quarterly                                 ###
### Last update: September 12, 2025                                   ###


### Packages 
library(lme4)              
library(lmtest)
library(sandwich)
library(stargazer)
library(xtable)
library(tidyverse)
library(dplyr)
library(stargazer)
library(ggplot2)
library(jtools)
library(ggeffects)
library(estimatr)
library(ggcorrplot)
library(ggpubr)
library(stringr)
library(RColorBrewer)
library(colorRamps)
library(mgcv)
library(miceadds)
library(texreg)
library(haven)
library(corrplot)
library(margins)
library(marginaleffects)
library(psych) 
library(RColorBrewer)
library(effects)
library(modelsummary)
library(flextable)

### Load data

df_repli <- read.csv("df_repli.csv")


### DESCRIPTIVES

# Fig 1: ## Self-identified Ideology by PID 

# Re-order PID cat from left to right
df_repli$pid_cat <- factor(df_repli$pid_cat, levels = c("NDP","Green", "Liberal", "Bloc Québécois", "Conservative", "PPC", "No party"))

p_ideo_pid <- ggplot(data=subset(df_repli, !is.na(SL_ideology_1) & !is.na(pid_cat)), aes(y = pid_cat, x = SL_ideology_1))+
  geom_boxplot()+
  labs(y= "Party Identification", x = "Ideological Self-Placement\n(Left to Right)") +
  theme_bw()

print(p_ideo_pid)
try(ggsave("Fig1.png", p_ideo_pid, width = 6, height = 5, units = "in", dpi = 600))


## Fig 2: Histogram of self-identified ideology  

df_repli$ideo_v2 <- as.factor(df_repli$SL_ideology_1)
levels(df_repli$ideo_v2 )[1] <- "0"
levels(df_repli$ideo_v2 )[11] <- "10"

total_obs <- nrow(df_repli[!is.na(df_repli$ideo_v2),])

d_ideo <- df_repli %>%
  filter(!is.na(ideo_v2)) %>%
  group_by(ideo_v2) %>%
  summarise(count = n()) %>%
  mutate(proportion = count / total_obs)

ideo_bar <- ggplot(d_ideo, aes(x = ideo_v2, y = proportion)) +
              geom_bar(stat = "identity", fill = brewer.pal(n = 11, name = "RdBu"), color = "black") +
              labs(title = "", x = "Self-identified Ideology\n(Left to Right)", y = "Percentage") +
              scale_y_continuous(labels = scales::percent_format(accuracy = 1), limits = c(0,0.3)) +
              theme_bw()

print(ideo_bar)
try(ggsave("Fig2.png", ideo_bar, width = 4, height = 4, units = "in", dpi = 600))


## Kernel Density Plot of policy-based ideology by Ideological self-identification

# bar plot of r-skew for centrist, leftists and rightists

d_mix5 <- df_repli %>%
  filter(!is.na(right_positions_p)  & !is.na(df_repli$ideo_orientation) ) %>%
  group_by(ideo_orientation, right_positions_p) %>%
  summarise(n = n()) %>%
  mutate(proportion = n / sum(n)) # proportions sum to 1 within ideological groups

d_mix5$right_positions_p_f <- as.factor(d_mix5$right_positions_p)

sum(d_mix5$proportion[d_mix5$ideo_orientation == "Left" & d_mix5$right_positions_p > 0.5]) # 11%
sum(d_mix5$proportion[d_mix5$ideo_orientation == "Right" & d_mix5$right_positions_p < 0.5]) # 43%
sum(d_mix5$proportion[d_mix5$ideo_orientation == "Centrist" & d_mix5$right_positions_p != 0.5]) # 75%
sum(d_mix5$proportion[d_mix5$ideo_orientation == "Centrist" & d_mix5$right_positions_p > 0.5]) # 33% of centrist lean on the right on issues
sum(d_mix5$proportion[d_mix5$ideo_orientation == "Centrist" & d_mix5$right_positions_p < 0.5]) # 42% of centrist lean on the left on issues

df_repli$ideo_orientation <- factor(df_repli$ideo_orientation,
                                    levels = c("Left", "Centrist", "Right"))


mix_density <- ggplot(df_repli[!is.na(df_repli$ideo_orientation),], aes( right_positions_p )) +
                geom_density(adjust = 3) +  
                facet_wrap(~ideo_orientation) +
                labs(x = "Policy-based Ideology\n(Left to Right)", y = 'Density') + # title: Right-leaning bias of centrist voters
                theme(axis.text.y = element_text(size = 12),  # Increase font size of y-axis tick labels
                      axis.text.x = element_text(size = 12),
                      axis.title.x = element_text(size = 13),
                      axis.title.y = element_text(size = 13),
                      strip.text = element_text(size = 13),
                      legend.title = element_text(size = 13),
                      legend.text = element_text(size = 12)) +
                theme_bw()

print(mix_density)
try(ggsave("Fig3.png", mix_density, width = 8, height = 4, units = "in", dpi = 600))

# N obs
nrow(df_repli[!is.na(df_repli$ideo_orientation) & !is.na(df_repli$right_positions_p),])

## Fig 4: Issue position by ideology

df_repli$SL_covid2[df_repli$SL_covid == "COVID-19 vaccines should not be mandatory for anyone"] <- "No one"
df_repli$SL_covid2[df_repli$SL_covid == "COVID-19 vaccines should be mandatory for vulnerable people only"] <- "Vulnerable"
df_repli$SL_covid2[df_repli$SL_covid == "COVID-19 vaccines should be mandatory for vulnerable people, travelers, and federal civil servants only"] <- "Vulnerable, travelers, civil servants"
df_repli$SL_covid2[df_repli$SL_covid == "COVID-19 vaccines should be mandatory for everyone"] <- "Everyone"

df_repli$SL_covid2 <- as.factor(df_repli$SL_covid2)


cols2 <- c("No one"  = brewer.pal(n = 4, name = "RdBu")[4],
           "Vulnerable" = brewer.pal(n = 4, name = "RdBu")[3],
           "Vulnerable, travelers, civil servants" = brewer.pal(n = 4, name = "RdBu")[2],
           "Everyone"  = brewer.pal(n = 4, name = "RdBu")[1])

covid_ideo <- ggplot(df_repli[!is.na(df_repli$SL_covid2) & !is.na(df_repli$ideo_orientation),], aes(x = ideo_orientation, fill = SL_covid2)) +
  geom_bar(position = "fill")+ # stacked bars add up to 1
  labs(title = "Covid", y = "Proportion that takes position", x = "Ideological self-placement", fill = "Position on\nmandatory vaccines") +
  scale_fill_manual(values = cols2,
                    labels = function(x) str_wrap(x, width = 20)) +
  guides(fill = guide_legend(reverse = TRUE,nrow = 4)) +
  theme(axis.text.y = element_text(size = 12),  # Increase font size of y-axis tick labels
        axis.text.x = element_text(size = 12),
        axis.title.x = element_text(size = 13),
        axis.title.y = element_text(size = 13),
        legend.title = element_text(size = 13),
        legend.text = element_text(size = 12)) 


df_repli$SL_climate2[df_repli$SL_climate == "Keep greenhouse gas emissions as they are"] <- "Keep the same"
df_repli$SL_climate2[df_repli$SL_climate == "Cut greenhouse gas emissions by 33% by 2050"] <- "Cut by 33%"
df_repli$SL_climate2[df_repli$SL_climate == "Cut greenhouse gas emissions by 66% by 2050"] <- "Cut by 66%"
df_repli$SL_climate2[df_repli$SL_climate ==  "Reach net-zero greenhouse gas emissions by 2050"] <-  "Net-zero"

df_repli$SL_climate2 <- as.factor(df_repli$SL_climate2)

cols2 <- c( "Keep the same" = brewer.pal(n = 4, name = "RdBu")[4],
            "Cut by 33%"= brewer.pal(n = 4, name = "RdBu")[3],
            "Cut by 66%" = brewer.pal(n = 4, name = "RdBu")[2],
            "Net-zero"= brewer.pal(n = 4, name = "RdBu")[1])

climate_ideo <- ggplot(df_repli[!is.na(df_repli$SL_climate2) & !is.na(df_repli$ideo_orientation),], aes(x = ideo_orientation, fill = SL_climate2)) +
  geom_bar(position = "fill")+ # stacked bars add up to 1
  labs(title = "Climate Change", y = "Proportion that takes position", x = "Ideological self-placement", fill = "Position on\nemissions") +
  scale_fill_manual(values = cols2,
                    labels = function(x) str_wrap(x, width = 20)) +
  guides(fill = guide_legend(reverse = TRUE,nrow = 4)) +
  theme(axis.text.y = element_text(size = 12),  # Increase font size of y-axis tick labels
        axis.text.x = element_text(size = 12),
        axis.title.x = element_text(size = 13),
        axis.title.y = element_text(size = 13),
        legend.title = element_text(size = 13),
        legend.text = element_text(size = 12)) 


df_repli$SL_immigration2[df_repli$SL_immigration == "Decrease immigration levels by 30%" ] <- "Decrease by 30%"
df_repli$SL_immigration2[df_repli$SL_immigration == "Decrease immigration levels by 10%"] <- "Decrease by 10%"
df_repli$SL_immigration2[df_repli$SL_immigration == "Increase immigration levels by 30%" ] <- "Increase by 30%"
df_repli$SL_immigration2[df_repli$SL_immigration == "Increase immigration levels by 10%"] <- "Increase by 10%"

df_repli$SL_immigration2 <- as.factor(df_repli$SL_immigration2)


cols2 <- c("Decrease by 30%"  = brewer.pal(n = 4, name = "RdBu")[4],
           "Decrease by 10%" = brewer.pal(n = 4, name = "RdBu")[3],
           "Increase by 10%" = brewer.pal(n = 4, name = "RdBu")[2],
           "Increase by 30%"  = brewer.pal(n = 4, name = "RdBu")[1])

immigration_ideo <- ggplot(df_repli[!is.na(df_repli$SL_immigration2) & !is.na(df_repli$ideo_orientation),], aes(x = ideo_orientation, fill = SL_immigration2)) +
  geom_bar(position = "fill")+ # stacked bars add up to 1
  labs(title = "Immigration", y = "Proportion that takes position", x = "Ideological self-placement", fill = "Position on\nimmigration") +
  scale_fill_manual(values = cols2,
                    labels = function(x) str_wrap(x, width = 20)) +
  guides(fill = guide_legend(reverse = TRUE,nrow = 4)) +
  theme(axis.text.y = element_text(size = 12),  # Increase font size of y-axis tick labels
        axis.text.x = element_text(size = 12),
        axis.title.x = element_text(size = 13),
        axis.title.y = element_text(size = 13),
        legend.title = element_text(size = 13),
        legend.text = element_text(size = 12)) 


df_repli$SL_deficit2[df_repli$SL_deficit == "Decrease deficit by 6% and cut back on social services"] <- "Decrease by 6%"
df_repli$SL_deficit2[df_repli$SL_deficit == "Decrease deficit by 2% and cut back on social services"] <- "Decrease by 2%"
df_repli$SL_deficit2[df_repli$SL_deficit == "Increase deficit by 2% to spend more on social services"] <- "Increase by 2%"
df_repli$SL_deficit2[df_repli$SL_deficit == "Increase deficit by 6% to spend more on social services"] <- "Increase by 6%"

df_repli$SL_deficit2 <- as.factor(df_repli$SL_deficit2)

cols2 <- c( "Decrease by 6%"   = brewer.pal(n = 4, name = "RdBu")[4],
            "Decrease by 2%"  = brewer.pal(n = 4, name = "RdBu")[3],
            "Increase by 2%"  = brewer.pal(n = 4, name = "RdBu")[2],
            "Increase by 6%" = brewer.pal(n = 4, name = "RdBu")[1])

deficit_ideo <- ggplot(df_repli[!is.na(df_repli$SL_deficit2) & !is.na(df_repli$ideo_orientation),], aes(x = ideo_orientation, fill = SL_deficit2)) +
  geom_bar(position = "fill")+ # stacked bars add up to 1
  labs(title = "Deficit", y = "Proportion that takes position", x = "Ideological self-placement", fill = "Position on\ndeficit") +
  scale_fill_manual(values = cols2,
                    labels = function(x) str_wrap(x, width = 20)) +
  guides(fill = guide_legend(reverse = TRUE,nrow = 4)) +
  theme(axis.text.y = element_text(size = 12),  # Increase font size of y-axis tick labels
        axis.text.x = element_text(size = 12),
        axis.title.x = element_text(size = 13),
        axis.title.y = element_text(size = 13),
        legend.title = element_text(size = 13),
        legend.text = element_text(size = 12))

# Combine plots

issue_ideo <- ggarrange(deficit_ideo, immigration_ideo, climate_ideo, covid_ideo, ncol = 2, nrow=2)
print(issue_ideo)

try(ggsave("Fig4.png", issue_ideo, width = 12, height = 9, units = "in", dpi = 600))

# N obs

nrow(df_repli[!is.na(df_repli$SL_covid2) & !is.na(df_repli$ideo_orientation),])
nrow(df_repli[!is.na(df_repli$SL_climate2) & !is.na(df_repli$ideo_orientation),])
nrow(df_repli[!is.na(df_repli$SL_immigration2) & !is.na(df_repli$ideo_orientation),])
nrow(df_repli[!is.na(df_repli$SL_deficit2) & !is.na(df_repli$ideo_orientation),])

### HYPOTHESIS TEST 

# H1

reg_h1_2 <- lm_robust(Y ~ ideo_prox*conjoint, data = df_repli, clusters = id) # not significant
summary(reg_h1_2)

modelsummary(
    reg_h1_2,
    stars = FALSE,     
    statistic = c("p.value"),    
    output = "regH1.xlsx")
  

## Fig 5: Plot the effect of ideo_dist on vote choice by conjoint condition
h1_pred_2  <- ggeffect(reg_h1_2 , terms = c("ideo_prox", "conjoint"))

p_h1_2 <- ggplot(h1_pred_2, aes(x = x, y = predicted,
                                color = group, linetype = group, fill = group)) +
          geom_line(size = 1.2) +
          geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
                      alpha = 0.2, color = NA) +
          labs(
            x = "Ideological proximity",
            y = "Probability of vote",
            color = "Experimental\ncondition",
            linetype = "Experimental\ncondition",
            fill = "Experimental\ncondition"
          ) +
          theme_bw() +
          theme(plot.title = element_blank()) +
          scale_color_manual(
            values = c(brewer.pal(n = 9, name = "Greys")[5], brewer.pal(n = 9, name = "Greys")[9]),
            labels = c("Policy info.\npresent", "Policy info.\nabsent")
          ) +
          scale_linetype_manual(
            values = c("solid", "dashed"),
            labels = c("Policy info.\npresent", "Policy info.\nabsent")
          ) +
          scale_fill_manual(
            values = c(brewer.pal(n = 9, name = "Greys")[5], brewer.pal(n = 9, name = "Greys")[9]),
            labels = c("Policy info.\npresent", "Policy info.\nabsent")
          )

  
print(p_h1_2)
try(ggsave("Fig5.png", p_h1_2, width = 6, height = 4, units = "in", dpi = 600))

## Table C2: include p-values and export as .tex

texreg(
  list(reg_h1_2),
  dcolumn = TRUE,
  booktabs = TRUE,
  include.ci = FALSE,
  stars = c(0),  # disable significance stars
  include.se = TRUE,   # show standard errors (in parentheses below coef)
  include.pvalues = TRUE,  # add separate p-value column
  use.packages = FALSE,
  label = "reg",
  caption = "Regression of Vote Choice and Candidate Rating on Ideological Proximity and Experimental Condition",
  float.pos = "H",
  file = "h1_reg.tex"  # exports to this .tex file
)


# Extract p-values from the regression model
p_values <- summary(reg_h1_2)$coefficients[, 4]

# Format p-values to 2 decimal places using sprintf
formatted_p_values <- setNames(sprintf("%.2f", p_values), names(p_values))

# Print the formatted p-values with variable names
print(formatted_p_values)


### H2

## Main model

reg_h2_2 <- lm_robust(Y ~  right_positions_p*ideo + ideo_prox, data = df_repli[df_repli$conjoint == "reduced",], clusters = id) 
summary(reg_h2_2)

modelsummary(
  reg_h2_2,
  stars = FALSE)


texreg(
  list(reg_h2_2),
  dcolumn = TRUE,
  booktabs = TRUE,
  include.ci = FALSE,
  stars = c(0),  # disable significance stars
  include.se = TRUE,   # show standard errors (in parentheses below coef)
  include.pvalues = TRUE,  # add separate p-value column
  use.packages = FALSE,
  label = "reg",
  caption = "Regression of Vote Choice on Ideological Proximity by Candidate Ideological Label (Reduced Conjoint)",
  float.pos = "H",
  file = "h2_reg.tex"  # exports to this .tex file
)


# Extract p-values from the regression model
p_values <- summary(reg_h2_2)$coefficients[, 4]

# Format p-values to 2 decimal places using sprintf
formatted_p_values <- setNames(sprintf("%.2f", p_values), names(p_values))

# Print the formatted p-values with variable names
print(formatted_p_values)


# Table 2: AME
m_h2_2  <- summary(margins(reg_h2_2, variables = c("right_positions_p"), at = list(ideo = levels(df_repli$ideo)))) # is slope for extreme left and extreme right candidates significantly different from 0? # yes

# Sample size
N <- reg_h2_2$nobs
N_c <- reg_h2_2$nclusters
R2 <- round(reg_h2_2$r.squared,2)
R2_adj <- round(reg_h2_2$adj.r.squared,2)

# Convert to xtable
tab <- xtable(m_h2_2,
              caption = "AME of Policy-based Ideology on Vote Choice by Candidate Ideological Label (reduced conjoint)" ,
              label = "ame_h2",
              align = c("l","l","c","c","c","c","c","c","c"))

# Add N as a row at the bottom
addtorow <- list()
addtorow$pos <- list(nrow(m_h2_2))

addtorow$command <- c(
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{N =", N, "}} \\\\ \n"),
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{N Clusters =", N_c, "}} \\\\ \n"),
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{R$^{2}$ =", R2, "}} \\\\ \n"),
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{Adj. R$^{2}$ =", R2_adj, "}} \\\\ \n")
)
addtorow$pos <- as.list(rep(nrow(tab), length(addtorow$command))) # Each command at the table end

print(tab, 
      include.rownames = FALSE,
      add.to.row = addtorow,
      sanitize.text.function = identity,
      caption.placement = "top")

# Figure 6
vote_min_pred_2 <-ggeffect(reg_h2_2, terms = c("right_positions_p", "ideo"))

cols3 <- c("10: Right" = blue2red(n = 5)[1],
           "7.5" = blue2red(n = 5)[2],
           "5" = "black",
           "2.5" = blue2red(n = 5)[4],
           "0: Left" = blue2red(n = 5)[5])

df_repli$ideo <- factor(df_repli$ideo,
                                    levels = c("0: Left", "2.5", "5", "7.5", "10: Right"))

p_vote_min_2 <- ggplot(vote_min_pred_2, aes(x, predicted, 
                                            colour = group, 
                                            linetype = group, 
                                            fill = group)) + 
  geom_line(size = 1.2) + 
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
              alpha = 0.1, colour = NA) +
  ylim(0.1, 0.7) +
  labs(x = "Voter Policy-based Ideology\n(Left to Right)", 
       y = "Probability of Voting for Candidate", 
       colour = "Candidate\nIdeological\nLabel",
       linetype = "Candidate\nIdeological\nLabel") +
  theme_bw() + 
  scale_colour_manual(values = cols3) +
  scale_fill_manual(values = cols3) +
  scale_linetype_manual(values = c("solid", "dashed", "dotdash", "twodash", "longdash")) +
  guides(fill = "none")

print(p_vote_min_2)
try(ggsave("Fig6.png", p_vote_min_2, width = 6, height = 5, units = "in", dpi = 600))


################################################################################
                        ### Supplementary Materials ###
################################################################################

## Conjoint attributes

# plot the effects of all attributes
cong_vote2 <- lm_robust(Y ~ cong_climate  + cong_covid + cong_deficit + cong_immig + ideo_prox + ethnic + gender, data = df_repli, clusters = id) # starting from male, observational data from pre-experiment survey
cong_vote3 <- lm_robust(Y ~ ideo_prox + ethnic + gender, data = df_repli[df_repli$conjoint == "reduced",], clusters = id) # rescale ideo var to allow reasonable comparison across attributes

# Extract coefficients and confidence intervals from cong_vote2
coefs_cong_vote2 <- tidy(cong_vote2)

# Extract the coefficient and confidence intervals for ideo_prox, gender, ethnicity from cong_vote3
ideo_prox_coef <- tidy(cong_vote3) %>%
  filter(term == "ideo_prox" | term == "ethnicWhite" | term == "genderWoman") %>%
  mutate(conjoint = "Reduced")  # Add a column for proximity type

# Filter the ideo_prox, gender, ethnicity terms from cong_vote2 and add proximity type
ideo_prox_cong_vote2 <- coefs_cong_vote2 %>%
  filter(term == "ideo_prox"|term == "ethnicWhite" | term == "genderWoman") %>%
  mutate(conjoint = "Full")

# Combine the two ideological proximity datasets
ideo_combined <- bind_rows(
  ideo_prox_cong_vote2,
  ideo_prox_coef
)

# Remove the conjoint column from coefs_cong_vote2 to avoid duplication when combining
coefs_cong_vote2 <- coefs_cong_vote2 %>%
  filter(term != "ideo_prox" & term != "(Intercept)" & term != "ethnicWhite" & term != "genderWoman")

# Combine all coefficients into one dataset
combined_coefs <- bind_rows(
  coefs_cong_vote2,
  ideo_combined
)

combined_coefs[1,10] <- "Full"  
combined_coefs[2,10] <- "Full"  
combined_coefs[3,10] <- "Full"  
combined_coefs[4,10] <- "Full"

# Update the term factor levels to maintain order
combined_coefs$term <- factor(combined_coefs$term, levels = c(
  "cong_climate",
  "cong_covid",
  "cong_deficit",
  "cong_immig",
  "ideo_prox",
  "ethnicWhite",
  "genderWoman"
))

# Add a new column for labeling proximity types
combined_coefs <- combined_coefs %>%
  mutate(label = case_when(
    term == "cong_climate" ~ "Climate\n(congruence)",
    term == "cong_covid" ~ "Covid\n(congruence)",
    term == "cong_deficit" ~ "Deficit\n(congruence)",
    term == "cong_immig" ~ "Immigration\n(congruence)",
    term == "ideo_prox" ~ "Ideological proximity",
    term == "ethnicWhite" ~ "Ethnicity\n(White vs. Chinese)",
    term == "genderWoman" ~ "Gender\n(Woman vs. Man)",
  ))

# Remove rows with NA in the conjoint column
combined_coefs <- combined_coefs %>%
  filter(!is.na(conjoint))

# Figure C1: Coefficients plot
p_cong_vote2 <- ggplot(combined_coefs, aes(x = estimate, y = term, color = conjoint, shape = conjoint)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), position = position_dodge(width = 0.5)) +
  labs(
    x = "Effect on Probability\nof Vote",
    y = "Attributes",
    color = "Conjoint",  # Legend title for color
    shape = "Conjoint"   # Legend title for shape
  ) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkgrey") + # Add a vertical line at x = 0
  scale_x_continuous(breaks = seq(-0.10, 0.15, by = 0.02)) +
  scale_y_discrete(labels = combined_coefs$label) +
  scale_color_manual(values = c("Full" = "black", "Reduced" = "blue")) + # Define colors
  scale_shape_manual(values = c("Full" = 16, "Reduced" = 17)) + # Define shapes: circle (16) and triangle (17)
  theme_bw()
  theme(
    axis.text.y = element_text(size = 12),  # Increase font size of y-axis tick labels
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(size = 13),
    axis.title.y = element_text(size = 13),
    legend.title = element_text(size = 14),  # Set legend title font size
    legend.text = element_text(size = 12),
    legend.position = "right"
  ) 

# Display the plot
p_cong_vote2

# Save the plot
try(ggsave("FigC1.png", p_cong_vote2, width = 8, height = 6, units = "in", dpi = 600))

# Table C1
texreg(
  list(cong_vote3, cong_vote2),
  dcolumn = TRUE,
  booktabs = TRUE,
  include.ci = FALSE,
  stars = c(0),  # disable significance stars
  include.se = TRUE,   # show standard errors (in parentheses below coef)
  include.pvalues = TRUE,  # add separate p-value column
  use.packages = FALSE,
  label = "reg",
  caption = "Regression of Vote Choice on Policy Congruence and Ideological Proximity",
  float.pos = "H"  # exports to this .tex file
)

# Extract p-values from the regression model (Full)
p_values <- summary(cong_vote2)$coefficients[, 4]

# Format p-values to 2 decimal places using sprintf
formatted_p_values <- setNames(sprintf("%.2f", p_values), names(p_values))

# Print the formatted p-values with variable names
print(formatted_p_values)

# Extract p-values from the regression model (Reduced)
p_values2 <- summary(cong_vote3)$coefficients[, 4]

# Format p-values to 2 decimal places using sprintf
formatted_p_values2 <- setNames(sprintf("%.2f", p_values2), names(p_values2))

# Print the formatted p-values with variable names
print(formatted_p_values2)

## For Table C2, see line 295 above

## Centrist Non-centrist Analysis

# convert centrist to factor
df_repli$centrist <- as.factor(df_repli$centrist)

h1_hetero_ideo4 <- lm_robust(Y ~ ideo_prox*conjoint*centrist, data = df_repli, clusters = id) # not significant
summary(h1_hetero_ideo4)
modelsummary(
  h1_hetero_ideo4,
  stars = FALSE)


h1_hetero_ideo_pred4 <- ggeffect(h1_hetero_ideo4, terms = c("ideo_prox", "conjoint", "centrist"))

filtered_data <- h1_hetero_ideo_pred4 %>%
  filter(!(facet == "Centrist" & x < -5))  # Assuming 'group' represents the panel and 'x' is the slope

levels(filtered_data$group)  <- c("Full", "Reduced") # Rename levels for clarity

#### Fig C2: plot the interaction
p_h1_hetero_ideo4 <- ggplot(filtered_data, aes(x = x, y = predicted,
                                               colour = group, 
                                               linetype = group, 
                                               fill = group)) +
                      geom_line(size = 1.2) +
                      geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
                                  alpha = 0.2, colour = NA) +
                      facet_wrap(~facet) +   # <<--- facet by your variable
                      labs(x = "Ideological Proximity", 
                           y = "Probability of Voting for Candidate", 
                           color = "Conjoint",
                           linetype = "Conjoint") +
                      scale_x_continuous(limits = c(-10, 0)) +
                      theme_bw() +
                      theme(
                        plot.title   = element_blank(),
                        axis.text.y  = element_text(size = 12),
                        axis.text.x  = element_text(size = 12),
                        axis.title.x = element_text(size = 13),
                        axis.title.y = element_text(size = 13),
                        strip.text   = element_text(size = 13),
                        legend.title = element_text(size = 13),
                        legend.text  = element_text(size = 12)
                      ) +
                      scale_color_manual(
                        values = c(brewer.pal(n = 9, name = "Greys")[5],
                                   brewer.pal(n = 9, name = "Greys")[9]),
                        labels = c("Policy info.\npresent", "Policy info.\nabsent")
                      ) +
                      scale_linetype_manual(
                        values = c("solid", "dashed"),
                        labels = c("Policy info.\npresent", "Policy info.\nabsent")
                      ) +
                      scale_fill_manual(
                        values = c(brewer.pal(n = 9, name = "Greys")[5],
                                   brewer.pal(n = 9, name = "Greys")[9]),
                        labels = c("Policy info.\npresent", "Policy info.\nabsent")
                      ) +
                      guides(fill = "none")   # hide redundant fill legend
                                   
  
print(p_h1_hetero_ideo4)
try(ggsave("FigC2.png", p_h1_hetero_ideo4, width = 8, height = 4, units = "in", dpi = 600))


# Table C3: Export table

texreg(
  h1_hetero_ideo4,
  dcolumn = TRUE,
  booktabs = TRUE,
  include.ci = FALSE,
  stars = c(0),  # disable significance stars
  include.se = TRUE,   # show standard errors (in parentheses below coef)
  include.pvalues = TRUE,  # add separate p-value column
  use.packages = FALSE,
  label = "reg",
  caption = "Regression of Vote Choice on Ideological Proximity, Experimental Condition and Centrist Self-Placement",
  float.pos = "H" # exports to this .tex file
)


# Extract p-values from the regression model
p_values <- summary(h1_hetero_ideo4)$coefficients[, 4]

# Format p-values to 2 decimal places using sprintf
formatted_p_values <- setNames(sprintf("%.2f", p_values), names(p_values))

# Print the formatted p-values with variable names
print(formatted_p_values)

# Table C4: Margins

m_h1_hetero_ideo4 <- summary(margins(h1_hetero_ideo4, variables = c("ideo_prox"), at = list(conjoint = levels(df_repli$conjoint), centrist = levels(df_repli$centrist)))) # is slope for extreme left and extreme right candidates significantly different from 0? # yes

# Testing diff in slope for centrists across conjoints (linear hypothesis tests)
coef_h1_hetero_ideo4 <-  coef(h1_hetero_ideo4)
coef_h1_hetero_ideo4[2]  # effect of ideo_prox for centrists in full conjoint
coef_h1_hetero_ideo4[2] + coef_h1_hetero_ideo4[5] # effect of ideo_prox for centrists in reduced conjoint
  ## Just need to check stat. sig of coef_h1_hetero_ideo4[5]: ideo_prox:conjointreduced 


# Sample size
N <- h1_hetero_ideo4$nobs
N_c <- h1_hetero_ideo4$nclusters
R2 <- round(h1_hetero_ideo4$r.squared,2)
R2_adj <- round(h1_hetero_ideo4$adj.r.squared,2)

# Convert to xtable
tab <- xtable(m_h1_hetero_ideo4,
              caption = "AME of Ideological Proximity on Vote Choice by Ideological Orientation and Conjoint" ,
              label = "ame_h1",
              align = c("l","l","c","c","c","c","c","c","c","c"))

# Add N as a row at the bottom
addtorow <- list()
addtorow$pos <- list(nrow(m_h1_hetero_ideo4))

addtorow$command <- c(
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{N =", N, "}} \\\\ \n"),
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{N Clusters =", N_c, "}} \\\\ \n"),
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{R$^{2}$ =", R2, "}} \\\\ \n"),
  paste("\\hline \n \\multicolumn{9}{l}{\\textit{Adj. R$^{2}$ =", R2_adj, "}} \\\\ \n")
)
addtorow$pos <- as.list(rep(nrow(tab), length(addtorow$command))) # Each command at the table end

print(tab, 
      include.rownames = FALSE,
      add.to.row = addtorow,
      sanitize.text.function = identity,
      caption.placement = "top")


### Post-hoc Power Analysis for H1

# --- 1. Define data structure ---
set.seed(123)

n_resp <- 983
n_obs <- 5898
id <- factor(rep(1:n_resp, length.out = n_obs))

# Simulate predictors
binary <- rbinom(n_obs, 1, 0.5)        # binary predictor
continuous <- rnorm(n_obs)             # continuous predictor
interaction <- binary * continuous

# Simulate outcome (true interaction effect = 0 for now)
y <- 0.1*binary + 0.1*continuous + 0.0*interaction + 
  rnorm(n_resp)[id] + rnorm(n_obs, sd = 1)

dat <- data.frame(id, binary, continuous, y)

# --- 2. Fit mixed model ---
m <- lmer(y ~ binary * continuous + (1|id), data = dat)

# --- 3. Power curve across effect sizes ---
effect_sizes <- seq(0, 0.20, by = 0.02)

results <- lapply(effect_sizes, function(es) {
  fixef(m)["binary:continuous"] <- es
  
  p <- powerSim(m, test = fixed("binary:continuous"), nsim = 100)
  ps <- summary(p)   # extract numeric results
  
  data.frame(effect_size = es,
             power = ps$mean,
             lower = ps$lower,
             upper = ps$upper)
})

pc_df <- do.call(rbind, results)

# export table
write.csv(pc_df, "power_curve_interaction2.csv", row.names = FALSE)

# Fig D1: Power Analysis Plot
power <- ggplot(pc_df, aes(x = effect_size, y = power)) +
        geom_line(size = 1.2, color = "darkblue") +
        geom_ribbon(aes(ymin = lower, ymax = upper),
                    alpha = 0.2, fill = "skyblue") +
        geom_hline(yintercept = 0.80, linetype = "dashed", color = "red") +
        scale_y_continuous(
          limits = c(0, 1),
          breaks = seq(0, 1, by = 0.2),        # major ticks every 0.2 (20%)
          minor_breaks = seq(0, 1, by = 0.1),  # minor ticks every 0.1 (10%)
          labels = scales::percent_format(accuracy = 1)
        ) +
      scale_x_continuous(
        limits = c(0, 0.2),
        breaks = seq(0, 0.2, by = 0.02),        
        minor_breaks = seq(0, 0.2, by = 0.01)
      ) +
        labs(x = "Interaction Effect Size (binary × continuous)",
             y = "Power",
             title = "Power Curve for Interaction Effect") +
        theme_bw() 


try(ggsave("D1.png", power, width = 6, height = 6, units = "in", dpi = 600))

################################################################################
                            ### End of File ###
################################################################################

